宏基因组实战2. 数据质控fastqc, Trimmomatic, MultiQC, khmer
本文英文原版见下方github链接，由中科院朱微金博士翻译、测试、并进行中文注释和补充，全网首发“宏基因组”公众号。
https://2017-cicese-metagenomics.readthedocs.io/en/latest/toc.html
前情提要
如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章
测试数据
刘博士帮助把测试数据建立了一个百度云同步共享文件夹，有非常多的好处，请读完下文再决定是否下载：
下载被墙的数据；很多数据存在google, amazon的部分服务器国内无法直接下载，而服务器一般科学上网不方便，下载数据困难。大家下载失败的数据请到共享目录中查找；
预下载好的软件、数据库；有很多需要下载安装、注册的软件(在线安装包除外)，其实已经在共享目录了，节约小伙伴申请、下载的时间；
数据同步更新；任何笔记或教程不可避免的有些错误、或不完善的地方，后期通过大家的测试反馈问题，我可以对教程进行改进。共享目录不建议全部下载或转存，因为文件体积非常大，而且还会更新。你转存的只是当前版本的一个备份，就不会再更新了。建议直接在链接中每次逐个下载需要的文件，也对文件有一个认识过程。
方便结果预览和跳过问题步骤；服务器Linux在不同平台和版本下，软件安装和兼容性问题还是很多的，而且用户的权限和经验也会导致某些步骤相关软件无法成功安装(有问题建议选google、再找管理员帮助；想在群里提问或联系作者务必阅读《如何优雅的提问》)。在百度云共享目录中，有每一步的运行结果，读者可以下载查看分析结果，并可基于此结果进一步分析。不要纠结于某一步无法通过，重点是了解整个流程的分析思路。
最后送上本教程使用到的所有文件同步共享文件夹链接：http://pan.baidu.com/s/1hsIjosk 密码：y0tb 。
数据质控
https://2017-cicese-metagenomics.readthedocs.io/en/latest/quality.html # 有时连接不稳定打不开，等会就会好。或访问它更早版本的链接如下：
https://2017-dibsi-metagenomics.readthedocs.io/en/latest/quality.html
安装软件
安装依赖关系
sudo apt-get -y update && \
sudo apt-get -y install trimmomatic python-pip \
samtools zlib1g-dev ncurses-dev python-dev unzip \
python3.5-dev python3.5-venv make \
libc6-dev g++ zlib1g-de安装 fastqc
wget -c http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip
unzip fastqc_v0.11.5.zip
cd FastQC
chmod +x fastqc
cd创建Python3.5虚拟环境
cd
python3.5 -m venv ~/py3
. ~/py3/bin/activate
pip install -U pip
pip install -U Cython
pip install -U jupyter jupyter_client ipython pandas matplotlib scipy scikit-learn khmer
pip install -U https://github.com/dib-lab/sourmash/archive/master.zip运行Jupyter Notebook
# 配置
jupyter notebook --generate-config -y
cat >>~/.jupyter/jupyter_notebook_config.py <<EOF
c = get_config()
c.NotebookApp.ip = '*'
c.NotebookApp.open_browser = False
c.NotebookApp.password = u'sha1:5d813e5d59a7:b4e430cf6dbd1aad04838c6e9cf684f4d76e245c'
c.NotebookApp.port = 8888
EOF
#
jupyter notebook &1. 测序数据准备
我们分析采用 Hu et al., 2016. 文章中数据的子集，下载数据
# 创建数据文件夹
mkdir data
cd data
# 下载测试数据
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_1.fastq.gz
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_2.fastq.gz
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249_1.fastq.gz
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249_2.fastq.gz
# 如果无法科学上网而下载失败，尝试在文提供的百度云中的data目录中下载
# 检查文件
md5sum *.gz
# 改原始文件为只读，防止被修改
chmod u-w *2. fastqc质量评估
# 质控所有gz压缩的原始数据，t启动多线程，一般与文件数量一致
fastqc *.gz -t 4
# 显示所有网页版质量评估报告文件，可下载本地或用firefox查看
ll *.html3. Trimmomatic去接头和低质量序列
下载Illumina双端接头序列
curl -O -L http://dib-training.ucdavis.edu.s3.amazonaws.com/mRNAseq-semi-2015-03-04/TruSeq2-PE.fa运行Trimmomatics
# 调用for循环批处理文件
for filename in *_1.fastq.gz
do
# 提取双端公共文件名，并输出检验
base=$(basename $filename _1.fastq.gz)
echo $base
# 运行去接头程序
TrimmomaticPE -threads 9 \
${base}_1.fastq.gz \
${base}_2.fastq.gz \
${base}_1.qc.fq.gz ${base}_s1_se \
${base}_2.qc.fq.gz ${base}_s2_se \
ILLUMINACLIP:TruSeq2-PE.fa:2:40:15 \
LEADING:2 TRAILING:2 \
SLIDINGWINDOW:4:2 \
MINLEN:25
done宏基因组拼接前必须去干净接头，防止引入人造序列对结果影响
4. 质控后再评估
fastqc *.qc.fq.gz -t 4
# 查看再次质控结果，与之前的比较试试
ll *.qc_fastqc.html图1. 比较质控前后第一个样品右端接头污染水平。上图质控前接头污染水平近10%，质控后接近0.
评估报告的结果非常多，自己多读读，不懂上fastqc官网看帮助。
5. MultiQC多样品报告汇总(可选)
需要python3.5
# 激活Pythone3环境
. ~/py3/bin/activate
# 安装包
pip install git+https://github.com/ewels/MultiQC.git
# 生成多样品报告
multiqc . #虽然是可选步骤，但对于多样品还是非常有意义的。可以方便比较，节省时间。
图2. 多样品质控前后比较。图像还是交互式的，鼠标悬停可显示样品名。
6. K-mer过滤
https://2017-cicese-metagenomics.readthedocs.io/en/latest/kmer_trimming.html
如果我们绘制样品k-mer丰度的柱状图，你会注意到存在大量的unqiue K-mers，即使测序质量很高，但它们也是由测序错误导致的。
图3. 序列末端低质量区有极高复杂度的kmer
本节继续在Python3下运行
# 对质控前后的数据统计单端丰度距离
abundance-dist-single.py -M 1e9 -k 21 SRR1976948_1.fastq.gz SRR1976948_1.fastq.gz.dist
abundance-dist-single.py -M 1e9 -k 21 SRR1976948_1.qc.fq.gz SRR1976948_1.qc.fq.gz.dist
# 只对高覆盖度中的低丰度kmer剪切(更可能是测序错误)；低覆盖度保留
interleave-reads.py SRR1976948_1.qc.fq.gz SRR1976948_2.qc.fq.gz | trim-low-abund.py -V -M 8e9 -C 3 -Z 10 - -o SRR1976948.trim.fq图4. kmer过滤原理：
只对高覆盖度中的低丰度kmer剪切(更可能是测序错误)；低覆盖度保
为什么要进行k-mer剪切
如果不做这步也是可以的。但会增加下游组装的工作量，本步可使结果更准确，并增加下游拼接速度，以及内存消耗。
unique-kmers.py SRR1976948_1.qc.fq.gz SRR1976948_2.qc.fq.gz
unique-kmers.py SRR1976948.trim.fq结果如下：
# 质控后的32-mers数据
Estimated number of unique 32-mers in SRR1976948_1.qc.fq.gz: 65344914
Estimated number of unique 32-mers in SRR1976948_2.qc.fq.gz: 85395776
Total estimated number of unique 32-mers: 112758982
# k-mer剪切后的数据
Estimated number of unique 32-mers in SRR1976948.trim.fq: 101285633
Total estimated number of unique 32-mers: 101285633结果只经过了简单的尾部过滤，k-mer的数量减少了10%以上，对下游分析的准确度和速度都非常有帮助。
按Kmer质控后的结果，感觉趣的再用fastqc评估一下，看看有什么变化？
接下来的文章来会介绍k-mer更大的用途，猜猜是什么？
猜你喜欢
热文：DNA提取发Nature
写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内五十位PI，六百多名一线科研人员加入。参与讨论，获得专业指导、问题解答，欢迎分享此文至朋友圈，并扫码加创始人好友带你入群，务必备注“姓名-单位-研究方向-职务”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决推荐生信技能树-微生物组版块(http://www.biotrainee.com/forum-88-1.html) 发贴，并转发链接入群，问题及解答方便检索，造福后人。
学习16S扩增子、宏基因组科研思路和分析技术，快关注“宏基因组”
点击阅读原文，跳转最新文章目录阅读